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Abstract 

We revisit and extend the original definition of discrete-time stochastic Petri nets, 
by allowing the firing times to have a “defective discrete phase distribution'*. We show 
that this formalism still corresponds to an underlying discrete-time Markov chain. The 
structure of the state for this process describes both the marking of the Petri net and 
the phase of the firing time for of each transition, resulting in a large state space. We 
then modify the well-known power method to perform a transient analysis even when 
the state space is infinite, subject to the condition that only a finite number of states 
can be reached in a finite amount of time. Since the memory requirements might still 
be excessive, we suggest a bounding technique based on truncation. 


! This research was partially supported by the National Aeronautics and Space Administration under 
NASA Contract No. NAS1-19480 while in residence at the Institute for Computer Applications in Science 
and Engineering (ICASE), NASA Langley Research Center, Hampton, VA 23681-0001. 
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1 Introduction 


In the past decade, stochastic Petri nets (SPNs) have received much attention from re- 
searchers in the performance and reliability arena and have been extensively applied to 
the performance and reliability modeling of computer, communication, manufacturing, and 
aerospace systems [4. 5, 7. 10, 23]. While there is agreement on the appropriateness of SPNs 
as a description formalism for a large class of systems, two radically different solution ap- 
proaches are commonly employed: simulation and state-space-based analysis. Simulation 
allows to associate general distributions to the duration of activities (SPN transitions), but 
it requires multiple runs to obtain meaningful statistics. This problem is particularly acute 
in reliability studies, where many runs might be required to obtain tight confidence intervals. 
With simulation, the state of the SPN is composed of the marking , describing the structural 
state of the SPN, and the remaining firing times, describing how long each transition in the 
SPN must still remain enabled before it can fire. The simulated time 0 is advanced by firing 
the transition with the smallest remaining firing time. 

State-space-based analysis has been mostly applied to SPNs whose underlying process is 
a continuous-time Markov chain (CTMC), that is, to SPNs with exponentially distributed 
firing times [3, 12, 25, 26]. Except for numerical truncation and roundoff, exact results are 
obtained, but the approach has two limitations: the number of SPN markings increases 
combinatorially, rendering unfeasible the solution of large models, and generally-distributed 
activities must be modeled using “phase-type (PH) expansion 7 ' [15]. PH distributions can 
approximate any distribution arbitrarily well, but it is difficult to exploit this fact in practice 
because the expansion exacerbates the state-space size problem. 

Discrete distributions for the timing of SPNs have received less attention. This is un- 
fortunate, since deterministic distributions (constants) are often needed to model low-level 
phenomena in both hardware and software, and the geometric distribution is the discrete 
equivalent of the exponential distribution and can approximate it arbitrarily well as the size 
of the step decreases. Furthermore, there is evidence supporting the use of deterministic 
instead of exponential distributions when modeling parallel programs [1]. 

If all the firing distributions are geometric with the same step, the underlying process is 
a discrete-time Markov chain (DTMC) [25]. Such SPNs can model synchronous behavior, 
as well as the main aspect of asynchronous systems: the uncertainty about the ordering 
of quasi-simultaneous events. A DTMC is described by a square one-step state transition 
probability matrix n and an initial state probability vector The state probability vector 
at step rt can be obtained with the iteration (power method ): 7rt n l = Tr^^n. This result 
was extended in [11] to include immediate transitions, which fire in zero time, and geometric 
firing distributions with steps multiple of a basic unit step, possibly with parameter equal 
one. that is, constants. [29] restates these results in more detail, and uses the concept of 
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weight to break the ties, following [3] and, more closely, [13]. Generalized Timed Petri Nets 
(GTPN) have also been proposed [19], where the steps of the geometric firing times for 
each transition can be arbitrary, unrelated, real numbers. A DTMC can be obtained by 
embedding, but the analysis is restricted to steady-state behavior and the state space of the 
DTMC can be infinite even when the underlying untimed PN has a finite reachability set. 
Analogous considerations hold for D-timed PNs [30]. 

We generalize and formalize the results in [13] and show how, using phase-expansion, a 
DTMC can be obtained even if the firing time distributions are not geometric, as long as 
firings can occur only at some multiple of a unit step. The state can then be described by 
the marking plus the phase of each transition. This extends the class of SPNs that can be 
solved analytically, but two limitations still exist: the existence of a basic step and the size of 
the state space. By using a fine step, arbitrary steps can be approximated, but this increases 
the state space. 

Approaches to solve models with a large state space have been proposed for both steady- 
state and transient analysis. [6] considers the reliability study of a SPN with exponentially 
distributed firing times, under the condition that the reachability graph is acyclic. The 
underlying CTMC is then acyclic as well, and a state can be discareded as soon as the 
transitions out of it have been explored, resulting in an algorithm offering large savings 
in memory and computations with respect to traditional numerical approaches. However, 
acyclic state spaces arise only in special cases, such as reliability models of non-repairable 
systems. 

For transient analysis of a general CTMC, Jensen’s method [21], also called uniformization 
[17, 27], is widely adopted. [18] outlines a dynamic implementation of the algorithm, where 
the state space is explored as the computation of the transient probability vector proceeds, 
not in advance, as normally done. This allows to obtain a transient solution even if the state 
space is infinite, provided that the transition rates have an upper bound. 

If the CTMC contains widely different rates, the number of matrix-vector multiplications 
required by uniformization can be excessive. Proposals to alleviate this problem are selective 
randomization [24] and adaptive uniformization [28], both based on the idea of allowing 
different uniformization rates, according to the set of states that can be reached at each 
step. The latter, in addition, can be used with infinite state spaces even if the rates have no 
upper bound. However, the method can incur a substantial overhead, and it appears that 
an adaptive step is advantageous only in special cases or for short time horizons. 
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In Sections 2, 3, and 4 we define the underlying untimed PN model, the class of DDP 
distributions used for the temporization of a PN, and the resulting DDP-SPN formalism, 
respectively. Section 5 discusses the numerical solution of a DDP-SPN, by building and 
solving its underlying stochastic process, a DTMC. Section 6, examines approaches to cope 
with large state spaces. 

2 The PN formalism 

We recall the (extended) PN formalism used in [12, 14]. A PN is a tuple (P, T. D ~ , D + , D ° , 
where: 


• P is a finite set of places, which can contain tokens. A marking p G IN^ defines 
the number of tokens in each place p G P, indicated by p p (when relevant, a marking 
should be considered a column vector). P~, D + . P°, and g are “marking-dependent”, 
that is, they are specified as functions of the marking. 

• T is a finite set of transitions. P fl T = 0. 

• Vp G P, Vf G T, V/t G IN |P| , D pt (p) G IN, D pt (p) G IN, and D° t {p) G IN are the 
multiplicities of the input arc from p to t, the output arc from t to p, and the inhibitor 
arc from p to t, when the marking is p, respectively. 

• >- C T x T is an acyclic (pre-selection) priority relation. 

• Vf G T, Vp G IN '' P \,g t (p) G {0.1} is the guard for t in marking p. 

• G IN^ P ' is the initial marking. 

Places and transitions are drawn as circles and rectangles, respectively. The number of 
tokens in a place is written inside the place itself (default is zero). Input and output arcs 
have an arrowhead on their destination, inhibitor arcs have a small circle. The multiplicity 
is written on the arc (default, is the constant 1); a missing arc indicates that, the multiplicity 
is the constant 0. The default value for guards is the constant 1. 

A transition t G T is enabled in marking p iff all the following conditions hold: 

1. g t {p) = 1. 

2. Vp G P, D~ t (p) < p p . 

3. Vp G P,D° tt (p) > p p or D° t {p) = 0. 

4. Vu G T, u t or u is not enabled in p (well defined because >- is acyclic). 
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Figure 1 

: Examples of DDP distributions. 


Let £{p) be the set of transitions enabled in marking p. A transition t £ £{p) can fire, 
causing a change to marking M(t, p), obtained from p by subtracting the “input bag'’ D~ t {p) 
and adding the “output bag v D, A (p) to it: M{t,p) = p — D, yt (p) + Dt t (p) — P + D,j(p), 
where D = D + -D~ is the incidence matrix. M can be extended to its reflexive and transitive 
closure by considering the marking reached from p after firing a sequence of transitions. The 
reachability set is given by Tl = {p : Ba € T* A p = M(a, p [0] )}, where T* indicates the set 
of transition sequences. 


3 Discrete phase distributions 

We now define the class B of (possibly defective) discrete phase (DDP) distributions, which 
will be used to specify the duration of a firing time in a SPN. A random variable X is said to 
have a DDP distribution, X ~ Z>, iff there exists an absorbing DTMC {A w : k £ IN} with 
finite state space A = {0. 1, . . .,«} and initial probability distribution given by [Pr{.4 [0] = 
?;}, i £ A], such that states A\ {0,n} are transient and states {0,n} are absorbing, and X is 
the time to reach state 0: X = min{A’ > 0 : A [k] - 0}. If Pr{yl [0] = 0} > 0, the distribution 
has a mass at the origin. If Pr{A [0] = i] > 0 and state i can reach state n, the distribution 
is (strictly) defective. 

V is the smallest class containing the distributions Const(O), Const(l), and Const(oo) 
and closed under: 

• Finite convolution: if X\ ~ and X 2 ~ T>, then X = X\ + X 2 ~ T>. 

• Finite weighted sum: if X\ ~ 2?, X 2 ~ 22 and B £ {0,1} is a Bernoulli random 
variable, then X = BX 1 + (1 — B) X2 ~ B. 

• Infinite geometric sum: if {X k ~ V \ k £ IN + } is a family of iid’s and N is a geometric 
random variable, then X = Yli<k<N X k ~ B. 
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Figure 2: Equivalent DTMC representations 

The geometric and modified geometric distributions with arbitrary positive integer step, 
Geom(a,u;) and ModGeom(a,u;), 0 < a < 1, u> £ IN + , the constant non-negative integer 
distribution, Const(u>), u £ IN, and any discrete distribution with finite non-negative integer 
support are special cases of DDP distributions. An example of a random variable with non- 
negative integer support which does not have a DDP distribution is N 2 , where N ~ Geom(o). 

Fig. 1 shows examples of DDP distributions. The “initial state” 6, for begin, has zero 
sojourn time and is introduced to represent graphically the initial probability distribution. 
We use this representation since it allows to estimate the “complexity” of a DTMC by 
counting the number of nodes and arcs in its graph. For simplicity, the last state, e.g., 
4 for Geom(a, 3) and 3 for Const(2), can be omitted if it is not reachable from b (if the 
distribution is actually not defective). Unfortunately, the DTMC corresponding to a given 
DDP distribution might not be unique, even if the number of states if fixed. For example, 
the time A' to reach state 0 for the DTMCs in Fig. 2, both with five nodes and seven arcs, 
has distribution Unif(0.3), that is. Pr{A' = ?'} = 1/4, for i £ {0, 1,2,3}. 

4 The DDP-SPN formalism 

SPNs are obtained when the time that must elapse between the instant a transition becomes 
enabled and the instant it can fire, or firing time, is a random variable. By restricting the 
firing times distributions to T>, we obtain the DDP-SPNs, corresponding to a stochastic 
process where the state has the form s = {g,<j>) £ IN |P ' x IN |T| . The structural component 
// is simply the current marking. The timing component 0 describes the current “phases”, 
the state for the DTMC chosen to encode the DDP distribution associated with the firing 
time of each transition. The firing time of a transition t elapses when its phase <p x reaches 0. 
Formally, a DDP-SPN is a tuple (P, T. D~, D + , D° , 4>, G, P, ^l°J, »-, ivj where: 

• (P, T, D~ , D + , D° , y, g, define a PN. 

• Vf £ T, Vp £ 71, $<(/t) C IN is the finite set of possible phases in which transition t can 
be when the marking is //. 


5 




• Vp e liyt 6 £(n),Vi,j e is the probability that the phase of t 

changes from i to j at the end of one step, when t is enabled in marking p. Hence, 

Z)je*tG) = I- 

• Vp € € £(p),Vf e T, Vi € $t(p),Vj € ^(jVt^p)),^,^, *,.;') is the prob- 

ability that the phase of t changes from i to j when u fires in marking p. Hence, 

• Vf 6 T, d>[°' 6 $t(p^) is the phase of t at time 0. 

• y- C T x T is an acyclic (post-selection) priority relation. 

• Vp € 11, VS C ^(p),Vf € 5,tt7 t |5(p) € IR + is the firing weight for t when S is the set of 

candidates to fire in marking p. 

A transition t € T is said to be a candidate (to fire) in state s = (p, <f>) iff all the following 
conditions hold: 

1. t £ £(p). 

2 . & = 0 . 

3. Vu € T, u f f or w is not a candidate in s (remember that »- is acyclic). 

Let C(s) be the set of candidates in state s. G«(p,», •) is the one-step transition probability 

matrix of the DTMC {<p^ : k £ IN}, with state space ^(aO, corresponding to the DDP- 
distributed firing time for transition t in marking p in isolation, that is, assuming that no 
other transition firing affects the firing time of t. However, if another transition u fires 
before f, leading to marking p', the phase </> t of t will change according to the distribution 
F u , t (p, <*,•). Furthermore, after the firing of ig the phase of t will evolve according to 
G t (p', •,•), which might differ from G ( (p,»,»), it can even have a different state space, 
$ t (p') instead of $<(p). 

We stress that pre-selection and post-selection have a different semantic. Only in the 
case of immediate transitions the two become equivalent. Assume that only t and u satisfy 
the input, inhibitor, and guard conditions in p. We have three options, resulting in three 
different behaviors: 

• Specify a pre-selection priority between them, for example t >- u, so that u will not be 
enabled when t is. This means that the phase 6 t of t evolves according to G ( (p,*,«), 
while <j> u does not. The same effect would be achieved using a guard g u (fi ) = 0. 
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• Specify no pre-selection priority, but a post-selection priority between them, for exam- 
ple t u. This means that the phases of both t and u evolve in fi. The first one to 
reach phase 0 will fire but, in case of a tie, t will be chosen. However, if 0 n = 0 when 
t fires and if jFi, u (^,0.0) = 1, u might be a candidate in the new marking, and fire 
immediately after t. 


• Specify neither a pre-selection nor a post-selection priority between them. Then, as in 
the previous case, t and u are in a race to reach phase 0, but a tie is now resolved by 
a probabilistic choice according to the the weights: w t \ {<.„}(//) and tc u |{t, u }(//), respec- 
tively, where w is a normalization of w to ensure that the weights of the candidates in 
a marking sum to one. 


Let be the state of the DDP-SPN at step n. Then, the process {(pM, d>I"I) : n G IN} 

is a DTMC with state space S C IN' P ' x IN |r| . Its one-step transition probability matrix 
n is determined by considering the possibility of simultaneous firings. Consider a state 
s = (//,d>). If C(s) ^ 0, one of the candidates will fire immediately, and the sojourn time in 
s is zero. Otherwise, the sojourn time in s is one. Following GSPN [3] terminology, we call 
s a vanishing or tangible state, respectively. Hence, s is tangible iff 6 > 0. 

Let S s . s ' be the set of possible event sequences events leading from a tangible state 
s = (//, 6 ) to a tangible state s' = (/d, <p ') in one time step: 


= [a = fi^K d> (1) , t^\ . . . /i (n-1 *, d> (n_1, f (n_1) . ^ (n) . 0 (n) ) : 

n > O.pl 0 * = = n' , <f>^ — 4 > , 

Vf € T(p.),G't(p,d> ( ,(ji>| 0) ) > 0, 

v ? ,o < i < tm (,) e C(n {i \d> {t) ),n {t+1] = M{t li \n li) ), 


(1) 

(2) 
(3) 


(1) considers the one-step evolution of the phases for the enabled transitions in isolation, 
while (2) and (3) consider the sequentialized firing in zero time of zero or more transitions 
at the end of the one-step period. Hence, (/d‘), d>hl) is a vanishing state, for 0 < ? < n. 

The value of the nonzero entries of n is obtained by summing the probability of all 
possible sequences leading from s to s': 


n. 


= e ( n G t {fi,4> t ,(f> { t o} ) 

<7€S ls , \te£M 

/ n — 1 


i=o \teT ) I 


In a practical implementation, n is computed one row at a time. The complexity of 
computing row s of n can be substantial, depending on the length and number of sequences 


7 




Figure 3: (0, 0) can reach an infinite number of markings in one time step. 

in U' S s , s >. If U S ' S SlS ' is infinite, special actions must be taken. This can happen for two 
reasons: 

• 71 is itself infinite, and state s can reach an infinite number of states in a single 
step. Consider, for example, a single queue with batch arrivals of size N > 0, where 
N ~ Geom(a), as in Fig. 3. Following the firing of f, a geometrically distributed 
number of tokens will be placed in p 2 : when the token is finally removed from pi (by the 
firing of v ), p 2 contains N tokens with probability a(l — a) ,,v_1 . This represents a batch 
arrival of size N at the server modeled by place p 2 and transition y. Unfortunately, 
finiteness of 71 is an undecidable question for the class of Petri nets we defined, since 
transition priorities alone make them Turing equivalent [2]. 

• S s . s ' can be infinite for a particular s'. If 71 is finite, this requires the presence of 
arbitrarily long paths over a finite set of vanishing states, just as for a “vanishing 
loop” in a GSPN [11]. In a practical implementation, these cycles can be detected and 
managed appropriately. 

The size of the DTMC underlying a DDP-SPN is affected by the choice of the representa- 
tion for the DDP distributions involved. Consider, for example, the DDP-SPN in Fig. 4(a), 
and assume that transitions G, t 2 , and f 3 have firing time distributions Const(l), Unif(0, 3), 
and Const(2), respectively. The corresponding DTMCs obtained using the two representa- 
tions of Fig. 2 for Unif(0,3) are shown in Fig. 4(b) and 4(c), respectively. The number of 
states is ten in the first case, seven in the second (the value of 4> t is specified as whenever 
t is not enabled and either it cannot become enabled again or its phase is going to be reset 
upon becoming enabled). The difference between the size of the two DTMCs is due to a 
lumping [22] of the states, and it would be even greater if t 3 had a more complex distribu- 
tion. By postponing the probabilistic decision as much as possible, the second DTMC lumps 
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Figure 4: The effect of equivalent Unif(0, 3) representations. 

states (011,«12), (Oil, #22). and (Oil, «32) of the first DTMC into a single one. (Oil, *32), 
and states (Oil, *11) and (Oil, *21) into (011,»21). 

5 Analysis of DDP-SPNs 

When using a SPN to model a system, a reward rate p M is associated to each marking p. 
Starting from {(p^, <f>W) : n € IN}- it is then possible to define two continuous-parameter 
processes: {y(9),9 > 0}, describing the instantaneous reward rate at time 9: y{9) = p M ie\. 
where p(9) = ^[«nax{n<01 ? an d {T'(0).0 > 0}, describing the reward accumulated up to time 
9, Y(9) = f° Plt{r) dT. 

We consider the computation of the expected value of y(9p) and Y(9p) for finite values 
of 9p . Let = Pr{s^ = s} be the state probability vector at time n. Once 

the state-space S corresponding to the initial state (/ff°l, 0^) has been generated, any initial 
probability vector over S can be used for the initial probability vector x^, there is no 
requirement to use a vector having a one in position (/<^, <£>t°l) and a zero elsewhere. From 
x^, we can obtain 7rl n l iteratively, performing n matrix-vector multiplications: 

X W = 7T [n ' 1] n (4) 

Since the DTMC can change state only at integer times, 7 r(0) = 7 r^ for 9 6 [n, n + 1). Practi- 
cal implementations assume that the state space is finite and that the transition probability 
matrix II is computed before starting the iterations. The following shows the pseudo-code 
to compute E[y(9p)\ and E\Y{9p)} with the “power method”: 

1. “compute S. II, and x^”; 

2. Y «- 0; 7T «- ttM; 

3. for n — 1 to [#fJ do 
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